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Ground-State Phase Diagram of Quantum Heisenberg Antiferromagnets 
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The S = 1/2 and 5 = 1 two-dimensional quantum Heisenberg antiferromagnets on the anisotropic 
dimerized square lattice are investigated by the quantum Monte Carlo method. By finite-size-scaling 
analyses on the correlation lengths, the ground-state phase diagram parametrized by strengths of 
the dimerization and of the spatial anisotropy is determined much more accurately than the previous 
works. It is confirmed that the quantum critical phenomena on the phase boundaries belong to the 
same universality class as that of the classical three-dimensional Heisenberg model. Furthermore, 
for S = 1, we show that all the spin-gapped phases, such as the Haldane and dimer phases, are 
adiabatically connected in the extended parameter space, though they are classified into different 
classes in terms of the string order parameter in the one-dimensional, i.e., the zero-interchain- 
coupling, case. 
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I. INTRODUCTION 

Recently many low-dimensional antiferromagnets with 
excitation modes separated from a ground state by 
a finite energy gap have been synthesized, and ef- 
fects of impurities or magnetic fields on those mate- 
rials have been investigated experimentally in relation 
to the impurity-induced long-range order (LRO) and 
magnetic-field-induced LRO. For example, there are the 
S — 1/2 quasi-one-dimensional (QID) Heisenberg antU 
ferromagnet (HAF) with bajid dimeri2iation, CuG(e03,El 
S = 1 QIdJiAF's, NENPJa NDMAzJa NDMApJfl and 
PbNi2V208,El and S* = 1 QID HAF's with bond alter- 
nation, NTEAPeI and NTENP.Q Those materials have 
attracted our interest since they reveal various aspects 
of the quantum phase transition between the quantum- 
disordered spin-gapped phase and the classical (Neel) 
long-range-ordered phase. 

Intrachain spin interaction with or without bond al- 
ternation and interchain interaction are considered to be 
the most basic ingredients to understand the quantum 
phase transitions mentioned above. More explicitly, they 
are expected to be modeled effectively by the quantum 
HAF on the anisotropic dimerized square lattice, which 
is described by the Hamiltonian: 

"H — 2_^ Szij' • S2i+lj- + Ct2_^ ^2i+lj' • ^2i+2,j (1) 



Here Si^j is the quantum spin operator at site (i,j) on 
the square lattice. The first two terms in r.h.s. repre- 
sent the one-dimensional antiferromagnetic (AF) Heisen- 
berg chains with alternating coupling constants, 1 and a 
(0 < a < 1), and the last term does the AF interchain 
exchange interaction (J' > 0). We choose the x-axis as 



«j 



being along the chain direction and the y-axis as in the 
perpendicular one. The bond arrangement of this model 
is shown in Fig. n^. 

The ground state of decoupled chains, i.e., J' = 0, has 
been |Well established. In particular, that of the S = 1/2 
chainB is the dimer state with a finite spin gap except 
for the uniform case (a = 1), which has a critical ground 
state. On the other hand, for the 5=1 chaia there exist 
two spin-gapped phases: the Haldane phaseO at a > etc 
and the dimer phase at a < ctc-tJ At the critical point 
a = ttc between these two phases the gap vanishp^j-The 
value of Uc has been estimated to be 0.5879(6)£jO In 
both cases, the critical point is_considered to belong to 
the Gaussian universality class .E2l 

For the AF LRO to appear the higher-dimensionality 
effect, i.e., the interchain interaction J', is indispensable. 
In most of the numerical works reported so far, the ef- 
fect of interchain coupling has been examined in certain 
approximated |Or perturbed ways. For example, Sakai 
and Takahashita estimated the critical strength, J^, for 
the uniform case {a — 1) by the exact diagonalization 
method for the intrachain interactions combined with the 
mean-field approximation for the interchain interaction, 

and obtained fj^^^^'"^^ = ^ ^^ '^'^^^"^^ - ^■'^25. More 
recently, Koga and Kawakamill3 investigated the S — 1 
model by the cluster-expansion method, and obtained 
,][. = 0.056(1) for a = 1. However, there have been only 
a very limited number of numerical works, in which both 
of the interchain and. .intrachain interactions are treated 
on an equal footing.liJ£3 Such numerical analyses are cer- 
tainly required, since the mean-field-like approximatioa. 
is not necessarily appropriate even in the QID regime.E3 
In the present paper, we report the results of quantum 
Monte Carlo (QMC) simulations-fai-using the continuous- 
imaginary-time loop algorithmliZrEj on the S = 1/2 and 
5 = 1 HAF model described by Eq. (|l). The present 
paper is organized as follows. In Sec. Hi the method 
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FIG. 1. Anisotropic dimerized square lattice with alternat- 
ing intrachain coupling of strength 1 (thick solid line) and a 
(solid line), and the interchain coupling of strength J' (dashed 
line) . 



of our numerical analyses is explained. In Sees. H and 
IV , the ground-state phase diagram parameterized by the 
strength of the bond alternation, a, and that of the in- 
terchain coupling, J', is determined precisely for S = 1/2 
and S" = 1, respectively. Especially, for the S" = 1 system 
with a = 1, we obtain J^ = 0.043648(8), which is consist, 
tent with 0.040(5) suggested by the recent QMC work,Ej 
but is much more accurate. Furthermore, both in the 
S — 1/2 and S — \ systems, the quantum phase transi- 
tions between the spin-gapped phases and the AF-LRO 
phase are confirmed to belong to the same universality 
class with that of the 3D classical Heisenberg model: the 
exponent of the correlation lengths \s v = 0.71(3) for 
S = 1/2 and i/ = 0.70(1) for 5 = 1, which coincideSriairly 
well with that of the latter model, v = 0.7048(30).B We 
also show the results on the correlation length and the 
gap in the spin-gapped phase. In Sec. 0, the topology 
of the phase diagram is discussed in detail based on the 
result of the present QMC calculation. We show that all 
the spin-gapped phases, such as the Haldane and dimer 
phases, are adiabatically connected with each other in 
the extended phase parameter space. This is in a sharp 
constrast to the strict ID case, in which the spin-gapped 
phases are classified into different-classes in terms of the 
so-called string-order parameter .E3 It is of interest that 
the ID spin-gap phases, which have different hidden sym- 
metry, are connected without encountering any singular- 
ity in the 2D phase diagram. The final section is devoted 
to the concluding remarks. 



II. METHOD 

We consider the system descrived by the Hamilto- 
nian (|l|) with 5=1/2 and S = 1. The real-space size is 
Lx X Ly and the inverse temperature, i.e., the imaginary- 
time size, is /3 = 1/T. Periodic boundary conditions 
are imposed in the x- and j/-directions. We use the 
continuous-imagiiaary-time loop algorithm with multi- 
cluster updateED'ES Especially, for the S — . 1. s y stem we 
adopt the subspin-representation technique ,113123 in which 



the S = 1 system is represented by an S = 1/2 system 
with special boundary conditions in the imaginary-time 
direction. By using these techniques, we can perform the 
simulatation up to L^ x Ly = 336 x 48 with /3 = 100 for 
the S* = 1 case without encountering any difficulty. 

The imaginary-time dynamical structure factor is de- 
fined by 



L.,Lyl3 ^ 



dtdt'e-''^-^" 



(2) 
^(*"*')(5f(t)5,^(i')), 



where q = (g^,, Qy) is the wave-number vector, Sf{t) is the 
^-component of the spin on site i at imaginary time r, 
and (• • •) is the thermal average. By using Sd{qx,<ly,'^), 
the staggered correlation length along the a;-directiai|iT-, 
is then evaluated by the second-moment methodjHj 



t _ ^ / 5'd(7r, 7r,0) 

^"^ ^ 2^ V Sd{TT + 2TT/Lx,7T,0) 



(3) 



The correlation length in the y-direction, ^y, and that in 
the imaginary-time (r) direction, ^t, which is related to 
the gap A by A = 1/£,t, are calculated similarly. Finally 
the staggered susceptibility, Xs, is evaluated by 



Xs 



Sd{n,n,0) 
1 'sr^ 



(4) 



LxLyP . . 



didt'e""-('"'~''^)(S'f(i)5'Kt')). 



All the structure fjaetors are calculated by using the im- 
proved estimators.Ea The period of 10^-10'^ Monte Carlo 
steps (MCS) is used for thermalization and that of lO'^- 
10^ MCS for the evaluation of physical quantities. 

Near the critical point (ac, ^c) of ^^^ ground-state 
transition, the correlation lengths diverge as 



r 



^r 



r^" = t- 



(5) 

(6) 



where t is the distance from the critical point and v is the 
critical exponent for the correlation length. Hjejje we have 
put z = 1 assuming the Lorenz invariance.E3 Furthata. 
more, the following finitc-size-scaling (FSS) formulaEZI 
holds near (ac, J^) and T = for systems with the fixed 
ratio Lx '■ Ly : (3, 

^x/Lx^f{tLl''^,LlT) = f{tLl/^), (7) 

and similar ones for S,y and S,t, and 

X. ^ L-l"^g{tLl'\LlT) = Ll''^ g{tLl''^). (8) 

Here / and g are scaling functions and 7 the exponent 
for Xs (~ t~'^). Note that LxT is put constant in the 
above equations. We assume a polynomial up to the 
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FIG. 2. The ct-J' phase diagram of the S — 1/2 system 
at zero temperature. The (spin-gapped) dimer phase and the 
AF-LRO phase are described by D and AF, respectively. The 
statistical error of each data point is smaller than the symbol 
size. The dotted and dashed lines denote J' = 1 and a = J' , 
respectively. 
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3. Finite-size scaling plot of the inverse gap, ^t, for 



FIG 
5 = 1/2 with J' 



1 and L^ : Ly : (3 = 1 : 1 : 1. The crit- 



ical coupling Oc and the exponent u are estimated as to be 
0.31407(5) and 0.71(1), respectively. The dashed line repre- 
sents the scaling function, which is approximated by a poly- 
nomial of order 2. 



second order for the scaling functions. By using least- 
squares fitting, we obtain the critical point (ac, J^) and 
the associated critical exponents v and 7. 

In addition, at some points in the spin-gapped phase, 
we explicitly evaluate the correlation lengths, ^^ and ^y, 
and the gap. A, at T = in the thermodynamic limit 
LxjLy — > 00. For this purpose we extrapolate the simu- 
lated data first to the ground state T — *■ and then to 
the thermodynamic limit Lx,Ly — > 00. 



III. RESULTS FOR S = 1/2 



particular, we obtain the critical point, ac — 0.31407(5), 
on the line J' = 1 (dotted line in Fig. 0), which is signifi- 
cantly smaller than their estimate etc — 0.398. More im- 
portantly, in the present simulation, the critical exponent 
ly on the transition points is evaluated as i^ = 0.71(1), 
which is consistent with-y ~ 0.7048(30) for the 3D clas- 
sical Heisenberg model.Ell The similar results ti4ve been 
obtained for the 2D 1/5-depleted HAF modelE3 On the 
other hand, KI concluded v = 1. The reason of the these 
discrepancies might be due to the smallness of the system 
sizes and the inverse temperature used in the study by 
KI. 



A. Ground-state phase diagram 

In Fig. H we show the ground-state phase diagram of 
the S — 1/2 system obtained by the FSS analysis ex- 
plained in the previous section. As an example of the FSS 
analysis, we show in Fig. Hthat of ^r against a for J' = 1 
(dotted line in Fig. H). The aspect ratio of the (2 + 1)- 
dimensional system is taken as, Lx '■ Ly : (3 = 1 : 1 : 1. 
By the least-squares fitting, the exponent v and the crit- 
ical coupling Uc are estimated as 0.71(1) and 0.31407(5), 
respectively. Here, the figure in parentheses denotes the 
statistical error (Icr) in the last digit. We also perform 
the same analyses on other lines in the a- J' plane, whose 
results are presented by the solid circles in Fig. 0. For 
example, on the line a — J' (dashed line in Fig. g), we 
obtain ttc = Jc = 0.52337(3) and v = 0.71(3). In the 
phase diagram we can see that the ground state of the 
chain (J' =i-flljs_the dimer state with a spin gap except 
for a = l,E3&E2l and that the region of the AF-LRO 
phase enlarges monotonically as J' increases. 

Our phase diagram is qualitatively the same as that 
of Katoh and Imada (KI),Efl but not quantitatively. In 



B. Correlation lengths and the gap 

We also evaluate explicitly the ground-state correlation 
lengths and the gap on some points in the dimer phase 
by using the dynamic structure factors. Unless the points 
are very close to the critical line, these quantities in each 
systems with L {— L^ = Ly) saturate to the ground-state 
values at temperatures we have simulated. For example, 
the T-dependences of these quantities are not to be dis- 
cernible at T = 0.05 and 0.01 for a = J' = 0.4 and 0.5, 
respectively. On the other hand, the L-dependence still 
remains in sizes we have calculated. The L-dependences 
of the ground-state spatial correlation lengths and the 
gap are shown in Fig. for a = J' — 0.5, which is close 
to the critical point ac = J'c = 0.52337(3). Their val- 
ues in the thermodynamic limit are estimated by fitting 
ik{L) to ^fe(L) = i -bexp{-cL), where k = x, y, or r, S. 
is the value in the thermodynamic limit, and b and c are 
fitting parameters. As a result, we obtain ^a: = 3.0089(9), 
^y = 2.2097(6) and A = 0.32261(4) for a = J' = 0.4 and 
^x = 11.998(9), S,y = 9.312(10), and A = 0.0913(2) for 
a = J' = 0.5. As a (= J') becomes smaller, i.e., the 
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FIG. 4. System-size dependence of the correlation lengths 
and the inverse gap at zero temperature for 5* — 1/2 with 
a = J' = 0.5. The extrapolated values are denoted by solid 
symbols. 



J' 



H-II 



AF 



D \ 



/' H-I 



(4) 



(3) 



^i_ 



ttc 1 N 



(1) 



a 



(2) 



FIG. 5. Schematic ground-state phase diagram of the S=l 
system obtained before the present QMC analysis. Filled 
circles are the points where the corresponding ground state 
is well-established: (1) the dimer phase (D), (2) a;-parallel 
Haldane phase (H-1), (3) the AF-LRO phase (AF), and (4) 
y-parallel Haldane phases (H-11). 



system becomes more distant from the critical point, A 
becomes larger and ^ smaller. 



IV. RESULTS FOR S = 1 

A. Overview on the phase diagram 

Before going into the QMC analysis, let us here sum- 
marize the ground-state phase diagram of the 5 = 1 sys- 
tem argued so far, which is shown in Fig. a. For some 
points in the phase diagram the ground state is well un- 
derstood by the previous theoretical and numerical stud- 
ies: (1) {a, J') = (0,0): The system consists of a set 
of the isolated antiferromagnetically-coupled spin pairs. 
The ground state is a trivial tensor product of dimer 



singlets sitting on each bond. (2) (a, J') = (1,0): The 
system consists of isolated x-parallel Haldane chains. (3) 
(a, J') = (l;!): The system is a uniform and isotroijiC| 
2D HAF. There exists an AF LRO in the ground state.El 
(4) J' ~ oo: The system consists of j/-parallel Haldane 
chains. Note that in this limit the value of a becomes 
irrelevant (see also discussions in Sec. M). 

In their analysis Jay the cluster expansion method, 
Koga and Kawakamillj derived three phases in the QID 
region, and they called the regions which includes point 
(1), (2) and (3) the dimer phase, the Haldane phase, and 
the AF-LRO phase, respectively. The region near point 
(4) is another Haldane phase. Therefore we call here the 
region which includes point (2) the Haldane I (H-I) phase 
and the one which includes the line (4) the Haldane II 
(H-II) phase. For the uniform systems with a = 1 the 
H-I and H-II phases are equivalent when we exchange the 
roles of the x-axis and the y-axis, and of J' and 1/ J'. 



B. Haldane-AF phase transition in the 
non-dimerized system 

To demonstrate our FSS analyses, let us begin with 
critical behavior near (1, J^), which separates the H-I and 
AF phases (J^ ~ 0.04 due to Ref. ^. We sweep J' near 
supposed J^ with a fixed to unity. In Fig. ^, ^x/Lx, 
S,y/Ly, and £,t/P of the systems with L^ = Ly ^ f3 = L 
(= 16, 24,- • •, 64) are plotted. As one sees immediately, 
the data suffer from quite large corrections to scaling, i.e., 
the crossing point of the scaled correlation lengths with 
two different L's clearly shifts to larger J' as the system 
size increases. We attribute these large corrections to 
the strong spatial anisotropy in the coupling constants 
(J' < 1). Indeed, the value of £,y at J' = 0.0435 is 
9.29(7) (see Fig. ^ (b)) even for the largest system size 
L = 64, which is quite smaller than those in the other 
directions {^x = 60.4(3) and Cr = 24.7(2)). Among the 
three correlation lengths, £,x is the largest: the growth 
of the correlation is dominated only by the system size 
in the a;-direction. This indicates that we need larger 
lattices, especially in the a;-direction, in order to perform 
a precise FSS analysis. 

To simulate larger lattices with minimal costs, we 
therefore optimize the aspect ratio, Lx '■ Ly : /?, as 
explained below. Expecting that the scaled correlation 
lengths, S,x/Lx, ^y/Ly, and ^t//3, become nearly equal 
with each other at the critical point, we set the as- 
pect ratio Lx : Ly : /? as 7 : 1 : 2 based on the data 
presented in Fig. H. With this ratio we simulate sys- 
tems with Lx = 168, 224, 280, and 336 and perform 
the FSS analyses. The raw data of i,x/Lx, S,y/Ly, and 
^r//3 with this aspect ratio are shown in Fig. 0. Now 
S,y — 12.02(5) even for the smallest system size {Lx = 168 
and Ly = 24) at J' = 0.0435. The ratios ^x/Lx, 
S,y/Ly and ^t//? in a common range of J' become nearly 
equal, and corrections to scaling become much smaller 
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FIG. 6. Plot of the correlation lengths, (a) £,x/Lx, (b) 
£,y/Ly, and (c) £,t/ P, as a function of J' for 5=1 with a = 1 
and Lx '■ Ly : P = 1 : 1 : 1. 
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FIG. 7. Plot of the correlation lengths, (a) (,x/Lx, (b) 
^yjLy, and (c) £,t/ P, as a function of J' for S = 1 with a = 1 



and Lx 



/3 = 7 : 1 : 2. 



than in the previous ones as we expected. The FSS 
plot for S,x is shown in Fig. ||. The resultant J[. and v 
are as follows: {J'^.v) = (0.043648(9), 0.69(1)) from ^^, 
(0.043649(8), 0.71(1)) from S,y, and (0.043648(7), 0.69(1)) 
from S,T- Averaging these three values we conclude with 



7 = 1.373(3). 



(11) 



and 



J^ = 0.043648(8) 



= 0.70(1). 



(9) 



(10) 



Fixing the value of J^ and v thus determined, we next 
perform the FSS analysis on the staggered susceptibility 
Xs- The raw data of Xs vs. J' are shown in Fig. O and 
the fitting result is shown in Fig. nfl. The latter yields 



The exponents v and 7 we have obtained for the S — 1 
system agaift^gree with those of the 3D classical Heisen- 
berg model.Eil 

The critical point obtained just above is consistent 
with the previous result by the methpd involving the 
mean-field approximation, J'^ > 0.025,t3 and also-with 
that by the recent QMC method, J'^ = 0.040(5) B On 
the other hand, the present result significantly differs 
from that obtained by the cluster-expansion method, 
J[. — 0.56(l)Jl3 The present FSS analyses on our exten- 
sive QMC results make us possible to obtain the critical 
value with much higher accuracy than the other methods. 



^ 




FIG. 8. Finite-size-scaling plot of the correlation length, 
Crc for 5 = 1 with Q = 1 and La; : Ly : /3 = 7 : 1 : 2. The 
critical coupling J^ and the exponent v are estimated as to 
be 0.043648(9) and 0.69(1), respectively, by the least-squares 
fitting. The dashed line represents the scaling function, which 
is approximated by a polynomial of order 2. 
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FIG. 9. Plot of the staggered susceptibility, Xs, as a func- 
tion of J' for S = 1 with Q = 1 and L^ '■ Ly : [3 = 7 : 1 : 2. 



C. Ground-state phase diagram 

In a similar way as described in the previous subsec- 
tion, we obtain other critical points on the a-J' phase 
diagram as shown in Fig. |ll|. First, the scaled correla- 
tion lengths, ^x/Lx, S.y/Ly, and ^r//3, are calculated up 
to Lx = 64 with either a or J' fixed. Sweeping J' or 
a with sufficiently high resolution, we regard a crossing 
point of these £,'s as the critical point. Note that the 
optimal aspect ratio depends strongly on the value of a 
and J'. However, we adopt Lx : Ly : f3 = I : I : 1 
for simplicity. Although the results thus obtained suffer 
from relatively larger systematic corrections than those 
presented in the last subsection, the absolute magnitude 
of the systematic error in the estimates should be still 
smaller enough than the symbol size in Fig. O. 

For some critical points the FSS analysis as in the pre- 
vious subsection is also carried out. We obtain the ex- 
ponents v and 7 which are consistent with Eqs. ( |lO| ) and 
(O), respectively. This supports that the quantum crit- 
ical phenomena in the S = I system also belong to the 
same universality class as that of the 3D classical Heisen- 
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FIG. 10. Finite-size-scaling plot of the staggered suscepti- 
bility, Xs, for 5=1 with a = 1 and L^ '■ Ly : (3 = 7 : 1 : 2. 
The values of J^ and u are fixed to 0.043648 and 0.70, respec- 
tively. The exponent 7 is estimated as to be 1.373(3) by the 
least-squares fitting. The dashed line represents the scaling 
function, which is approximated by a polynomial of order 2. 



berg model. An exception is-the ID critical point located 
at (cKc, Jc) — (0.5879(6), 0),llil which separates the dimer 
phase from the Haldane phase. The apparent value of h' 
starts to deviate from ( p^ ) when a becomes closer to ac- 
This is attributed to the crossover to the critical phie-i 
nomena belonging to the Gaussian universality class .1121 
We confirm that the AF-LRO phase exists between the 
two spin-gap phases at least down to J' = 0.01 at a = ac. 
Although in the present simulation it is quite difficult to 
prove the existence of the AF-LRO phase at smaller J', 
we believe that the the point (ac, 0) is tricritical: the ID 
critical point is unstable against an infinitesimal inter- 
chain coupling and the AF LRO immediately amijears as 
the same as in the S = 1/2 uniform chain.ElErEl 

Interestingly, the H-II and D phases are adiabatically 
connected with each other. The gapless AF-LRO phase 
does not touch the line of a = as seen in Fig. p2, 
where the part of the whole phase diagram (Fig. ^l|) near 
a = is magnified. Indeed, on the a = line, which 
corresponds to the S* = 1 two-leg ladder, it is shown that 
there exists no critical point by the recent QMC study.cl 
Thus, the D phase can be identified with the H-II phase, 
and there are only two distinct spin-gap phases, H-I and 
H-II in Fig. |ll|. The closeness of the critical line to the 
a = line is due to the strong AF fluctiiations, which 
already exist in the two-leg ladder system.o 



D. Correlation lengths and the gap 



We obtain the explicit values of £,x, £,y, and A in the 



ground state at (a, J') = 
for systems with sizes Lx 



(1,0.04). They are calculated 
= 168, 224, 280, and 336 and 



with the aspect ratio L^^ : Lj, = 7 : 1 at T regarded 
as zero temperature. Their T-dependences are negligible 
at T = 0.01. We extrapolate the finite-size data to the 
thermodynamic limit in the same way as explained for 




FIG. 11. Ground-state phase diagram of the 5 = 1 system. 
The J' ajcis is logarithmically scaled for J' > 0.1 for conve- 
nience. The statistical error of each data point is smaller than 
the symbol size. 
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FIG. 13. Schematic ground-state phase diagram in the ex- 
tended parameter space for (a) S = 1/2 and (b) S = 1. The 
shaded areas denotes the AF-LRO phase on the ABC and 
BCD planes. The cross section on the plane R = 0.5 (x-y 
isotropic plane) is drawn based on the numerical results. 
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FIG. 12. Ground-state phase diagram of 5* = 1 system in 
the strongly-dimerized region (a < 0.05). 



the S = 1/2 system. We obtain ^^ = 39.2(1), ^y = 
5.67(1), and A = 0.0632(2). As J' becomes smaller, ^^ 
becomes smaller and A larger to reach at J' = tlie. 
single chain values ^^ = 6.0153(3) and A = 0.41048(6). E3 



V. DISCUSSIONS 

The analyses presented in the preceding section re- 
vealed that the ground-state phase diagram of the S — 1 
system has a rather complicated topology, i.e., the H- 
II and D phases are adiabatically connected with each 
other, though the channel between them is quite nar- 
row (Fig. ^l]). On the other hand, as for the H-I and 
D phases, in the ID system (J' = 0) these two spin- 
gapped phases are distinctivelMpSeparated by the criti- 
cal point at ac = 0.5879(6),BEI andUthey are distin- 
guished by the string-order parameter ,Ej which is zero in 
the former phase and finite in the latter one. The tran- 
sition can be viewed as a rearrangement of dimcr-singlet 



pattern betwften-4he (1,1)- and (2,0)-valence-bond-solid 
(VBS) states.E3H We emphasize that once J' is intro- 
duced, however, the string-order parameter should van- 
ish even- in the H-I phase, being similar to the 5 = 1 
ladder.a Still one may consider that the two phases es- 
sentially differ with each other since they are separated 
by the AF-LRO phase. If, however, we introduce the 
bond alternation also in the y-direction, the two phases 
can be connected without passing the gapless state as 
explained below. 

Let us consider the 2D HAF model defined in the ex- 
tended parameter space: 



Ti- — Jx \ 2_^ ^2ij" • S2i+1 J -I- ttx 2_^ ^2 



■j+l,j ■ '^2i+2,j 



«J 



/ . S«,2j 
i,3 



^,3 



5i,2j + l 



'Es 



i,2j + l ■ ^i.2j+2 



«,J 



(12) 



The original Hamiltonian (1^) corresponds to the case 



Jy 



J', and ay 



1. The 



with Jx — I, ctx 
Hamiltonian (03) is invariant under the exchange be- 
tween {Jx,ax) and {Jy,ay). 

To draw the phase diagram in this extended param- 
eter space, it is convenient to introduce a parameter, 
R = Jy/{Jx + Jy)- The limits Jy ^ and Jy ^ oo 
correspond to i? = and 1, respectively. Since ax (ay) 
becomes irrelevant in the limit J^, — > (Jx —>■ c)o), the 
whole phase diagram in the three-dimensional parame- 
ter space is shaped as a tetrahedron. In Fig. |3| we 
present the ground-state phase diagram parametrized by 
R, ax{l - R), and ayR for S = 1/2 (a) and S* = 1 (b). 
It should be noted that the phase diagram should be in- 
variant under the transformation, (i?, ax{^ — R), ayR) ^-> 



(1 — R, ayR, ax{l — R)), reflecting the symmetry in the 
Hamiltonian explained above. In the phase diagram, 
the edge AB (CD) corresponds to isolated x-parallel (y- 
parallel) decoupled chains, the edge AD isolated four-spin 
plaquettes, and the edge AC (BD) the two-leg ladders in 
y-direction (a;-direction). 

The face ABC (and also CDB) in Fig. |l^ corresponds 
to the original phase diagram shown in Figs, g and |ll|, 
though the y-parallel-chain limit J' — > oo in the origi- 
nal diagram is represented by one vertex C (B) in the 
new ones. In the extended phase diagram the shaded 
(unshaded) area represents the AF-LRO (spin-gapped) 
phase on the ABC- and CDB-faces. 

It should be emphasized that on the ACD- and DBA- 
faces there is no AF-LRO phase, since the system is one- 
dimensional dimerized two-leg ladder. There exist only 
the ID critical points discussed already. Especially, there 
is no critical point on the edge AD. Therefore, in the S — 
1 case, the three spin-gapped phases, H-I, D, and H-II, 
are connected by the path C-^A^D^B. Similarly, in the 
S — 1/2 case, the two dimer phases, which correspond to 
the vertex A and D, respectively, are connected directly 
by the path A^D. Thus, in both cases, there are only two 
phases, namely, the spin-gapped phase and the AF-LRO 



VI. CONCLUDING REMARKS 

In this paper, we have investigated the ground-state 
phase diagram of 5 = 1/2 and 5 = 1 HAF on the 
anisotropic dimerized square lattice by means of the ex- 
tensive QMC simulation with the continuous-imaginary- 
time loop algorithm and the FSS analyses. It is con- 
firmed that, for both S =1/2 and 1, the quantum critical 
phenomena in the model belong to the same universality 
class as that of the 3D classical Heisenberg model, ex- 
cept for the ID critical points, which belong to the same 
universality class as that of the Gaussian model. We 
have also demonstrated that the spin-gapped phases of 
the ID chain are connected when we introduce the inter- 
chain couplings with bond alternation. In the 2D system, 
only one spin-gapped phase exists in both of the S = 1/2 
and S — 1 systems. 

The results obtained in the present work are consid- 
ered to be the proper basis for investigation of peculiar 
phenomena observed in the QID HAF's mentioned at the 
beginning of this paper. We have already reported the 
QMC analysis on the site-dilution-induced AE LRO in 
these materials based on the Hamiltonian (|l|) £3 In order 
to discuss various experimental results quantitatively, it 
is certainly necessary to take into account other ingredi- 
ents than in Eq. (nl) , such as the next-nearest-neighbor in- 
trachain interaction and the single-ion anisotropy. They 
are beyond the scope of the present work. The present 
results, however, demonstrate the role of the higher di- 
mensionality which has been overlooked so far. 



Most of the numerical calculations in the present work 
have been performed on the DEC Alpha, SGI ORI- 
GIN 2000, SGI 2800, and RANDOM at the Materials 
Design and Characterization Laboratory, Institute for 
Solid State Physics, University of Tokyo and on the Hi- 
tachi SR-2201 at the Supercomputer Center, University 
of Tokyo. The program used in the present simulation 
was based on the library 'Looper version 2' developed by 
S.T. and K. Kato and also on the 'PARAPACK version 
2' by S.T. The present work was supported by the "Re- 
search for the Future Program" (JSPS-RFTF97P01103) 
of the Japan Society for the Promotion of Science. S.T's 
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